// 版权所有2009 Go作者。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

// 原始C代码和下面的长注释是
// 来自FreeBSD的/usr/src/lib/msun/src/e_sqrt.C和
// 随此通知而来。go代码是原始C.
// 
// ================================================================================================
// 版权所有（C）1993由Sun Microsystems，Inc.保留所有权利。
// 
// 由SunPro开发，SunPro是Sun Microsystems，Inc.的一家企业。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。
// /===============================================================
// /
// /\uIEEE754\uSQRT（x）
// 返回正确四舍五入的sqrt。
// ---------------------------------------
// |如果您有硬件sqrt，请使用硬件sqrt | 
// ---------------------------------------
// 方法：
// 使用整数算术的逐位方法。（速度慢，但携带方便）
// 1。标准化
// 在[1,4]中用2的偶数幂将x缩放到y:
// 找到一个整数k，使1<=（y=x*2**（2k））<4，然后
// sqrt（x）=2**k*sqrt（y）
// 让q=sqrt（y）在二进制点后截断为i位（q=1，
// i 0 
// i+1 2 
// s=2*q，y=2*（y-q）。（1）为了从q中计算q，我们要检查q是否为40
// 
// （i+1）2 
// （q+2）<=y.（2）
// i 
// -（i+1）
// /If（2）如果为false，则q=q；否则q=q+2。
// i+1i+1i 
// 
// /通过一些代数运算，不难看出
// （2）等价于
// （i+1）第52页
// i i 
// 以下递推公式计算：
// 如果（3）为假
// 
// s=s，y=y；（4）
// i+1 i+1 i 
// 
// 否则，
// /-i-（i+1）
// /s=s+2，y=y-s-2（5）
// i+1 i+1 i 
// 
// 人们可以很容易地使用归纳法来证明（4）和（5）。注意，（3）的左侧只包含i+2位，因此
// 不必做完整的（53位）比较
// in（3）.
// 3.最终取整
// 生成53位结果后，我们再计算一位。
// 与余数一起，我们可以确定
// 结果是否准确，大于1/2ulp还是小于1/2ulp 
// （它永远不会等于1/2ulp）.
// 可以通过检查
// 大+小是否等于大，以及大-小是否为
// 对于某些浮点数“大”和“小”是否等于大来检测舍入模式.
// 
// 
// 注释：舍入模式检测省略。常量“mask”、“shift”、
// 和“bias”在src/math/bits中找到。go 

// Sqrt返回x的平方根。
// 
// 特殊情况为：
// Sqrt（+Inf）=+Inf 
// Sqrt（±0）=±0 
// Sqrt（x<0）=NaN 
// Sqrt（NaN）=NaN 
func Sqrt(x float64) float64 {
	if haveArchSqrt {
		return archSqrt(x)
	}
	return sqrt(x)
}

// 注意：在某些系统上，Sqrt是在汇编中实现的。
// 其他程序集存根跳转到下面的func sqrt。
// 在Sqrt是一条指令的系统上，编译器
// 可以将直接调用转换为直接使用该指令。

func sqrt(x float64) float64 {
	// 特殊情况
	switch {
	case x == 0 || IsNaN(x) || IsInf(x, 1):
		return x
	case x < 0:
		return NaN()
	}
	ix := Float64bits(x)
	// 规范化x 
	exp := int((ix >> shift) & mask)
	if exp == 0 { // 次规范化x 
		for ix&(1<<shift) == 0 {
			ix <<= 1
			exp--
		}
		exp++
	}
	exp -= bias // 无偏指数
	ix &^= mask << shift
	ix |= 1 << shift
	if exp&1 == 1 { // 奇数exp，双x为偶数
		ix <<= 1
	}
	exp >>= 1 // exp=exp/2，平方根指数
	// 逐位生成sqrt（x）
	ix <<= 1
	var q, s uint64               // q=sqrt（x）
	r := uint64(1 << (shift + 1)) // r=将位从MSB移动到LSB 
	for r != 0 {
		t := s + r
		if t <= ix {
			s = t + r
			ix -= t
			q += r
		}
		ix <<= 1
		r >>= 1
	}
	// 最终舍入
	if ix != 0 { // 余数，结果不精确
		q += q & 1 // 根据额外位
	}
	ix = q>>1 + uint64(exp-1+bias)<<shift // 有效位+偏置指数
	return Float64frombits(ix)
}
